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ABSTRACT 

Elliptic equations in exterior regions frequently require a boundary 
condition at infinity to ensure the well-posedness of the problem. Examples 
of practical applications include the Helmholtz equation and Laplace’s 
equation. Computational procedures based on a direct discretization of the 
elliptic problem require the replacement of the condition at infinity by a 
boundary condition on a finite artificial surface. Direct imposition of the 
condition at infinity along the finite boundary results in large errors. A 

r 

sequence of boundary conditions is developed which provides increasingly 
accurate approximations to the problem in the infinite domain. Estimates 
of the error due to the finite boundary are obtained for several cases. 
Computations are presented which demonstrate the increased accuracy that can 
be obtained by the use of the higher order boundary conditions. The examples 
are based on a finite element formulation but finite difference methods can 
also be used. 
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I. 


Introduction 


Elliptic problems in exterior regions arise in many branches of physics. 
For example, the flow of an incompressible irrotational fluid about a 
body is described by the Laplace equation (e.g. Lamb [19]) 

A u = 0. 

The same equation arises in the study of electrostatics exterior to given 
surfaces (e.g. Stratton [30]). A different example is the exterior scattering 
problem for either acoustics or electromagnetism. In this case one wishes to 
solve the Helmholtz equation 


A u + k 2 u = 0 , 

with either Dirichlet or Neumann data specified on the bodies (Bowman, et al. 
[6], Muller [22]). For inhomogeneous media, k is a given function of the 
position. For some applications in plasma physics, k can be a nonlinear func- 
tion. 

In these cases, infinity can be regarded as a separate boundary. A 
condition at infinity is required to make the exterior problem well-posed. 

For the Laplace equation it is sufficient to impose a condition of regularity 
at infinity. In three dimensions this is 

(1.1) U = 0(i) r -*■ » , 

where r is the distance from a fixed (but arbitrary) origin (Kellogg [15]). 
For the Helmholtz equation one can impose the Sommerfeld radiation condition 



(Sommerfeld [30]) 


(1.2) - i k u = o(-^) 

A more exact form of (1.2) is 

(1.3) lim 

r -> oo 

y| = r 

where the integral is over spherical shells centered at r = 0 (Rellich 
[27], Rellwig [13]). The radiation condition ((1.2) - (1.3)) states that 
the solution corresponds to outgoing waves (see Wilcox [35], [36], 
for more details) . 

In many instances one is interested in problems with variable coefficients 
that approach a constant state at infinity. An extension of the theory of 
radiation conditions to problems with variable coefficients was developed by 
Vainberg [33]. The techniques to be described are valid for the variable 
coefficient case provided the coefficients approach constants at infinity 
at a sufficiently rapid rate. 

A numerical solution of an elliptic problem in an exterior region must 
be able to incorporate the radiation condition at infinity within the compu- 
tational procedure. Solution techniques based on eigenfunction expansions 
or asymptotic methods automatically accomplish this by the proper choice of 
expansion functions. When the free space Greenes function which satisfies 
the radiation condition is known, (e.g. constant coefficients) the difficulty 
of imposing the radiation condition can be avoided by reformulating the 
problem as a Fredholm integral equation. Such formulations for the Helmholtz 
equation can be found in Chertock [9], Kleinman and Roach [16], and Burton 
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and Miller [8]. Poggio and Miller [24] discuss a vector integral equation 
for the reduced Maxwell equations. Bayliss [4] shows that one can greatly 
increase the efficiency of the integral equation formulation by introducing 
an appropriate coordinate transform along the body. Schneck [28] has 
developed codes to solve the resultant integral equations for a range of ap- 
plications. 

The integral equation approach has several deficiencies. These methods 
are generally restricted to the constant coefficient case and require the 
inversion of a full matrix. This can result in storage difficulties, especially 
for three dimensional problems or problems with high frequencies. In addition, 
for many applications the matrix elements are expensive to compute [4]. Many 
mesh points are required to resolve the singularity in the kernel even for 
low frequencies. 

For the Helmholtz equation, an additional difficulty with the integral 
equation formulation is the possibility of interior resonances. It is well 
known that for certain values of k the integral equation becomes singular 
([9], [16]). These resonances are connected with eigenvalues of associated 
interior problems. Various attempts have been made to overcome this diffi- 
culty, but they generally increase the complexity of the integral equation 
approach (Ursell [32]). It will be shown that the proper formulation of 
radiation conditions can eliminate the possibility of eigenvalues. 

An alternative to the integral equation method is to couple an interior 
solution with a global functional of the solution on an artificial boundary. 

The global functional can be obtained by integral formulas using the free 
space Green's function or by using an expansion, typically obtained from 
separation of variables, to represent the solution exterior to the artificial 
boundary. Marin [21], Zienkiewicz, et al. [37] and others have studied 


- 3 - 



methods based on an integral relation over the artificial boundary, while 
Fix and Marin [10] have used a boundary condition based on separation of 
variables to solve problems in underwater acoustics. 

These methods incorporate the exact radiation condition at the cost of 
a non-local boundary condition. A disadvantage of these methods is that the 
non-local coupling over the artificial boundary is equivalent to the full 
matrix that would be obtained from the integral equation. Furthermore, for 
the Helmholtz equation, spurious eigenvalues can also occur with this formu- 
lation. Goldstein [11] has suggested extending the boundary conditions of 
Engquist and Majda to the elliptic case. However, this method is restricted 
to only a range of frequencies depending on the expansion parameter. No 
calculations using this method have been carried out to date. 

The method to be presented develops a sequence of local boundary condi- 
tions that are extensions of (1.1) and (1.2). These boundary conditions are 
then applied at a finite artificial boundary. As the order of accuracy of 
the boundary operator increases, the order of the highest derivative appearing 
in the boundary operator will also increase. 

The artificial surface will generally be assumed, for the proofs, to be 
the sphere r = r^. The resulting elliptic problem is then discretized and 
solved in the bounded region between the body and the artificial surface. The 
proposed boundary conditions are asymptotic in 1/r. Hence, for a given 
accuracy one can bring the artificial boundary further in when using the 
higher order boundary conditions. 

In many applications the solution is required only in the vicinity of 
the body. The far field solution can be calculated by a quadrature formula, 
such as Green T s formula, once u and t— are known along the body. For 
problems in potential flow one frequently is only interested in the solution 
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on the body. Because of this, we shall stress the accuracy of the boundary 

conditions as r^ decreases. Furthermore, the error estimates will be for 
2 

surface L errors. These errors are due to the imposition of the generalized 
radiation conditions at a finite boundary. 

These boundary conditions are related to a family of boundary conditions 
developed by Bayliss and Turkel [5] for time dependent problems. The boundary 
operators are differential relations which match the solution to an expansion 
in 1/r which is valid in a neighborhood of infinity. The boundary conditions 
generally involve derivatives of order greater than or equal to the order of 
the differential equation. Hence, these boundary conditions are different 
from the usual boundary conditions encountered in elliptic theory. 

For the Laplace equation the resultant discretization can be solved by 
fast iterative methods leading to a substanial improvement over the integral 
equation methods. For Helmholtz type equations, the error will have a depen- 
dence on the number of wave lengths between the body and the artificial 
surface. In several test cases it has been possible to constrict the compu- 
tational region as k increases. All numerical results in this study were 
obtained with a finite element program based on a band Gaussian solver. This 
has large storage requirements which limited the investigation of high fre- 
quencies. This storage can be reduced by the use of iterative 
methods. An iterative method for the Helmholtz equation based on a decay law 
for a corresponding hyperbolic problem was developed by Kriegsman and Morawetz 
[17]. Kriegsman and Morawetz have also implemented time dependent boundary 
conditions similar to those proposed here (private communication). Brandt [7] 
and Nicolaides [23] have developed iterative methods based on the multi-grid 
algorithm. These methods require Gaussian elimination on a coarse grid which 
provides some resolution of the solution. Hence, these methods have limited 
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use for high frequencies. The scheme of Nicolaides can be combined with 
the finite element method described in the appendix. Extensions of the 
fast solvers based on capacitance methods are also feasible (Proskorowski 
and Widlund [25]). 

In section 2 we develop the extension of the standard radiation con- 
ditions for both the Helmholtz and Laplace equations. Error estimates for 
Laplace* s equation are given in section 3 and for the Helmholtz equation in 
section 4. These chapters can be skipped by those only interested in the 
computational procedure. In section 5 computational results are presented 
for both the Laplace equation and the Helmholtz equation. For the Helmholtz 
equation both constant and variable indices of refraction are considered. 

The details of the finite element procedure are given in the appendix. 

II. Construction of Radiation Boundary Conditions 

Let u be an outgoing solution to the three dimensional Helmholtz 
equation 

(2.1) A u + k^u = 0 , 


exterior to a sphere r = r^. It is known (Atkinson [2], Wilcox [34]) 
that u can be represented by a convergent expansion 


( 2 . 2 ) 


u 


ikr 


e 


kr 


f 

j=0 (kr) j 


Here 0,(J> denote the angular variables of an (r,0,(j)) spherical coordinate 
system. The series (2.2) is uniformly and absolutely convergent and can be 
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differentiated term by term any number of times (Wilcox [34]). Fj (0 » ) > 
j > i can be obtained from the radiation pattern, Fq» formula 


Fj (9,4)} 


1 

(2i) j j! 


j 

n [iia-i) + Q]F n (e,(j)) 


9 


where Q is the Beltrami operator in the angular coordinates 9 and 4> 


- 1 9 , . fl 3 x . 1 9_ 

Q " sin9 39 ^ Sin0 30 ) + . 2 „ .,2 


sin 0 3(j) 


The Sommerfeld radiation condition for any such solution is 


-i k u + u^ = o(— ) 


(r -*■ °°) 


In fact it is clear from the expansion (2.2) that 


(2.3) 


-i k u + u 


= oo\) 

r 


(r . 


In numerical computations for which the exterior region is truncated at 
a finite value of r, say r = r^, a possible boundary condition to impose is 


(2. A) 


-i k u + u 


= 0 


r = r. 


However, this condition is very inaccurate and in fact it can easily be seen 
that it is not exact even for the first terra in the expansion (2.2). 

We will develop here a sequence of linear differential operators 
which provide more accurate extensions of the condition (2.4) by annihilating 
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the first m terms in the expansion (2.2). Thus the condition that at 
r = r 1 the solution lies in the null space of the operator B m can be 
considered as a procedure to match the solution to the first m terms in 
the expansion (2.2). 

An example of such an operator is 

B l u = ( l7 _ik+ r )u ' 


It is easily verified that 


ikr 


kr 


F(M) 


- 0 




for any function F (6 ,(J>) . Furthermore, it can be verified that for any 
function u having the expansion (2.2) 


(2.5) 



To develop more accurate conditions we consider the sequence of operators 


(2.6) B m - n (fj-ik+^1 = (fe-ik+^Vi 


3=1 


A straightforward calculation verifies that 


V 0 ■ 


for any function p of the form 
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«,= 1 


_ e ikr F(9,(!)) 
P “ 0 

(kr)* 


y • • * 9 


m 


Furthermore, if u is any function having an expansion of the form (2.2) 
then 


(2.7) 


B u 
m 


= h = 0(- 


2m+l 


r = r- 


It is clear that the leading order term in h will involve only the 

term of order m+ 1 in the expansion (2.2) and thus the boundary condition 

B u = 0 will match the solution to the first m terms in the expansion 
m 

(2.2). For any fixed k the errors in the boundary condition B^u will 
decrease at a faster rate (in (r^) 1 ) as r^-*-®. This is thus analogous to 
the use of higher order difference approximations where the errors decrease 
at a faster rate in a mesh size h. Thus one can expect a significant in- 
crease in efficiency by applying the higher order boundary conditions. We 
point that the individual terms in (2.2) are not solutions to (2.1) (unless 
k = 0) however the series (2.2) provides a description of the behavior of 
the solution as r + °° which is generally not the case for expansions based 
on complete sets of solutions to (2.1) (see Aziz and Kellogg [3]). We fur- 
ther point out that the operators B^ in (2.6) can be written as 


( 2 . 8 ) 


(-ik + 


8 vm 

9r 


+ lower order terms 


and thus they can be considered as generalizations of the Sommerf eld radia- 
tion condition. 

In the study of the Helmholtz equation we are interested in the be- 
havior of the errors in different parameter ranges. Specifically, 
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1. The error for fixed k and m as the position of the artificial 
surface, i.e. r^, varies. 

2. The error for fixed k and r^ as the order of the boundary 
operator m increases. 

3. The error for fixed r^ and m as the frequency k increases. 

Case 2 is of greater importance than case 1 as it is generally more 

expensive to increase the size of the computational domain than to implement 
the higher order boundary conditions. Results to be established in sections 
3 and 4 demonstrate that the higher order boundary condition can lead to 
significantly increased accuracy when applied at a fixed r^. 

The constants involved in the order relations (2.5) and (2.7) will, 
in general, depend on k. This is because for arbitrary problems the angu- 
lar functions F^(0,<f)) can be expected to grow with k. Thus as k increases 
the errors in (2.5) and (2.7) will not, in general, be bounded uniformly in 
k. However, the quantity k r^ is a natural non-dimensional quantity, which 
is in fact just the number of wave lengths to the artificial boundary r = r^. 
The error obtained by taking only a fixed number of terms in the series (2.2) 
can be expected to depend significantly on this quantity. A consequence of 
this is that the radial resolution will not increase, or will increase slowly 
as k increases. Numerical results illustrating this will be presented in 
section 5. 

We also consider the exterior Laplace equation. In this case it is 
well known that the analog of (2.2) holds, i.e., solutions have the multi- 
pole expansion 

i ^ V 9 ’^ 

(2.9) u = - £ J -j ~ • 

j=0 r 
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It then follows that the differential operators obtained by setting k = 0 
in (2.6) 


( 2 . 10 ) 


B - n (f + a=i) 

m 3r r 


> 


exactly annihilate the first m terms in (2.9). Similarly for any u 
with the expansion (2.9) 


( 2 . 11 ) 


B u 
m 


r = r„ 


o(- 


— ± — ) 
2m+r 
r 


Both (2.6) and (2.11) are valid for equations with variable coefficients 

provided the coefficients approach a constant state sufficiently rapidly 

at infinity. As indicated in the introduction we will be mainly interested 

in computing the solution on some surface. The solution at far fields points 

can then be obtained by a quadrature based on Green's formula. In the next 

two sections we prove theorems that demostrate that in some cases the expected 
2 

surface L errors are achieved. 

For computational ease one can replace all radial derivatives beyond 
the first by tangential (angular) derivatives. This is done by using the 
differential equation and is especially important for finite element applica- 
tions. The exact forms of B^ and B 2 are given by (5.2). 

The previous discussion has concentrated on the three dimensional Helm- 
holtz and Laplace equations. For the two dimensional Helmholtz equation, the 
solution u has the convergent expansion (Karp [14]) 


( 2 . 12 ) 


u-H 0 «cr) 

j=o r j=0 


where Hq and H^ are the Hankel functions of the first kind of order 0 


and 1. 


- 11 - 



As this expansion is difficult to work with we use the asymptotic expansion 

t 6 ] 


(2.13) 


nr i(kr-f)^,£ <e) 

n '/*r e 

j=0 


The boundary conditions based on (2.13) are 


m 


(2.14) 


3 

b = n (|- + 

m . - 9r r 

3 = 1 


- ik ) 


which are analogous to (2.6). 

We have concentrated on the homogenous Laplace and Helmholtz equations. 
The same boundary conditions can be used for the inhomogenous equation 


(2.15) 


A u + k u = F 


The estimates obtained in the next two sections apply to (2.15) provided 
that F decays sufficiently rapidly for large r. 

Boundary conditions based on the operators are nonstandard because 

of the high order of the derivatives involved. The regularity of the solutions 
up to the boundaries, using B^u =0 at the artificial surface, is guaranteed 
by the Agmon, Dougalis, Nirenberg theory [1]. Using the results of Lopatinskii 
[20] we need only consider the half plane problem without any lower order 
terms. Thus, we have reduced the problem to 



Au = 0 

x _> 0 

- 00 < y, z < 00 , 

(2.16) 

9 u = 0 

x = 0 



3x n 
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The only solution to (2.16) that decays as x 00 is u = 0. Hence, 
the complimentary condition of [1] is satisifed and regularity follows 
(see also [31] ) . 

When estimates for the error are obtained it is important that the 
constants in the estimation do not grow too rapidly as the artificial 
boundary approaches infinity. It is also important, especially for the 
Helmholtz equation, to eliminate the possibility of spurious eigenvalues 
for the resulting interior problem. Indeed we wish to obtain estimates 
that are uniform in k and do not have poles at discrete values of k. 

Such estimates will be obtained in the next two sections. 

III. Error Bounds for Laplace's Equation 

2 

We wish to obtain estimates for the L surface error that occurs 
when the regularity condition at infinity is replaced by the condition 
B^u =0 at a finite boundary. In this section we concentrate on Laplace's 
equation. As discussed in the introduction the main interest is in the 
errors that occur in the solution and its derivative along the inner boun- 
dary. 

For concreteness we shall consider the Neumann problem. All results 
are equally valid for the Dirichlet problem. 



Figure 1. 



The problem that we wish to solve is 


(3.1) 


A u = 0 

in 

3u 

8n g 

on r. 

u = 0(b 

f “*■ 00 


We replace this by the problem 


in Q 
on r 
on 

(see figure 1) . 

Let w be the error, w = u - v. Then w satisfies 


A v = 0 



B v = 0 
m 


(3.3) 


A w = 0 in 


9w n r 

= 0 on r i 


B w = B u = h on Y 

m m 


If r 2 is the sphere r = r^ then by (2.6) h ° ( 2m+F* 

r i 

In this section we will only consider the case that the artificial 
surface is the sphere r = r^. We will consider two cases. In 

theorem (3.1) we discuss the case that the body T ^ is also a sphere. In 
theorems (3.2) and (3.3) we consider general bodies, but only treat the 
boundary conditions B^ and B 2 * This is not a major restriction since B^ 
and B 2 are the boundary conditions of the greatest practical importance. 
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Since we are interested in surface errors we introduce the notation 


( 3 . 4 ) H w ll(r) = ff l w (y)| 2<1A y • 



Introducing spherical coordinates (r,0,<j>) (3.4) becomes 


w 


2 

(r) 


2t r 7T 


■a 


w(r,0,<j>) | ^r^sin 0 d0d(J> 


For simplicity we assume axial symmetry so that w is independent of (j). 
All the results are independent of this assumption. 

For the first part we assume that the body is also a sphere. The 
coordinates are scaled so that the surface is the sphere r = 1. We then 
have 


Theorem (3.1) 

Let g(0) be smooth and hence satisfy 


(3.5) 



(0) sin 0 d 0 


< co. 


Let w be the solution to (3.3). Then there exists a constant C, indepen- 
dent of m, such that 


(3.6a) 


(3.6b) 


| w 



n m+1 
Cr 

2m+l 

-i 


i 



II 


(r) 


Cr 
2m+l 
1 


1 < r < ^ 
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Note: 


(1) The bounds on the derivatives are necessary for the Dirichlet problem. 

(2) It follows from (3.6a) that 


W "(r ) ~ 0( -^> 
1 r. 


W "(l) ° ( 2nrKL ) 


Hence, the smallest errors occur along the inner surface. This is analogous 
to the Saint-Venant ! s principle. Lax [19] has considered decay laws for 
volume norms for positive elliptic operators. 

(3) Since C is independent of m the estimate shows that we can improve 
the accuracy by fixing r^ and increasing the order of the boundary operator. 


Proof 


We first consider the case that Dirichlet data is imposed, i.e., (3.1b) 
is replaced by u = g on By assumption (3.5) and using axial symmetry 

we have that 

oo 

(3.7) g (0) = ^ d . P. (cos 0 ) , 

j=0 J J 


where P . 

J 


are the Legendre polynomials. 


It then follows that 



g 2 (0) | sin 0 


d0 


oo 


E d 

j=0 


2 

j 


If g(0) is sufficiently smooth it is also known that for any n (Gottlieb 
and Orszag [12]) 
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c 

n 


(j+l) n 


where C n depends on n and g but not on j. Using separation of 
variables, the solution u to (3.1) is given by 


(3.8) 


d. < 

-T 


(3.9a) 


(3.9b) 


OO 

E d P (cos 0 ) 
j=0 


J+l 


fr (r ’ 0 ) - 


-E 

j=0 


(j+l)d^Pj (cos 0 ) 

J+2 


It is easy to verify, by induction, that the boundary operators given 

by (2.10) satisfy 


(3.10) 


B (r k) = 
m 


m 

n 

£=1 


(k+&) r k_m = A, r k-m 
k,m 


Clearly A, =0 for -m < k < 0. 
k,m — 

It follows from (3.2b), (3.9) and (3.10) that 


(3.11) 


B (w) = B (u) 


OO 


-E 


d .P. (cos 0 )A . 

-J_J _=J >3 

r j+m+l 


By separation of variables the error w has the form 


(3.12) 


w(r,0) = ^2 qjl-j Ceos 0)[r j - -j^} 

j=o r 
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The operator is applied to (3.10) and the result is evaluated 

at r = r^. Using the completeness of the Legendre polynomials, this 
series can be equated term by term, with (3.11). We then have 


(3.13) 


! 0 3 0 ,..., ni” 1 

44^ - 4 j - 4 1 

J 


j > m 


From the definition of A . it is easily verified that 

-J 


-(j+1) ,m 


J 


> 1 


and so 


(3.14) 


and 


I q j I i 


^ 1 


j > m , 




j+i 


•> i 




i!fi! r 


m 


r 2j 1 - 1 rj m 1 - 1 


1 ^ r <_ 
j > m 


Since, ||w|| 2 (r) = | q (r j - -~-) | 2 r 2 , the result (3.6a) follows, 

j r 

By differentiating the series for w (3.12), and repeating the process, (3.6b) 
follows . 

The proof for the case of Neumann data is similar. In this case (3.9) 
and (3.12) are replaced by 


d p (cos 0) 

u(r ’ e) ■ - E ^ 


j=0 


(j+D r j+l 


and 
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»(r,0) - 2 Vj <co,0)tr]+ ] • 

J-0 

We next consider the case where (3.1) is to be solved exterior to an 
arbitrary domain We restrict the proofs by only considering m = 1,2. 

The proof relies on the generalized maximum principle discussed by Protter 
and Weinberger [26]. To apply these methods we assume that has the 

property that every x € lies on the boundary of a ball contained in 
the exterior of F^. We first need 

in ft , 

on r 1 , 

on F 2 for j = 1,2 

Suppose there exists functions z^, z 2 such that 

A z^ <_ 0 in ft , 

3z i < 

(3.16) 0 on , 

BjZ 1 >_ h on T 2 j = 1,2 , 

and 
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in ft , 

on r , 

b j z 2 1 h on r 2 j = 1,2 , 

then 

£ u £ in ft 

We make several observations 

1. The normal derivative is taken in the direction away from the origin. 

3u 9u 

Hence, for a sphere 

2. Similar results hold when Dirichlet data is imposed on the inner body. 

3. The results are valid for any uniformly elliptic equation in ft. 

Proof . 

For j = 1 the lemma is a restatement of theorem 12 of [26]. When 
is used the result follows from the observation that the proof of 
Protter and Weinberger only requires that when w has a positive maximum 
on P 2 then B 2 w is positive. This follows from the form of B 2 . 

Based on this lemma we have 

Theorem 3.2 

For an arbitrary body V ^ and using the boundary operator B^ on 
the error satisfies 

| | . c pointwise. 

(3.18) M < ~2 
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Furthermore, using the operator B 2 on V 2 the error satisfies 


(3.19) | w | pointwise 


Proof . 

The boundary condition B^w can be expressed as 



c 

We choose for the lemma z^ = — and z 2 = -z., . Straightforward algebra 
shows that z^ and satisfy (3.16) and (3.17) respectively. (3.18) 

then follows from the conlcusion of the lemma. For the boundary condition 

C 

B, ; wo choose z^ = — ^ and z 2 = -z^ and again (3.19) fallows. 

r i 

We note that the exact solution to the exterior problem decays at 
least as fast as 1/r. Hence, the uniform error bounds given by (3.18) and 
(3.19) express a smaller relative error at the body r then near the arti- 
ficial surface r = r^. For some applications one wishes error bounds for 
the normal derivative on The bounds can be obtained by deriving a 

r\ 

Fredholm equation of the second kind for 7 ^ on T^. A theorem valid for 
k £ 0 will be given in the next section. 
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IV. 


Error Bounds for the Helmholtz Equation 


The problem that we wish to solve is 


(4.1a) 


(4.1b) 


(4.1c) 


2 

A u + k u = 0 exterior to 

3u 

3n ® on » 

9u . , ,1. 

3 r “ * * u _ °( r ) as r -*■ °o . 


As before we replace condition (4.1c) by 


(4. Id) 


B v = 0 
m 


on r 2 , 


where v also satisfies (4.1a, 4.1b) . Let w be the error, w=u-v 
Then w satisfies 


(4.2a) 

(4.2b) 

(4.2c) 


Aw + 


7 2 
k w 


3w 

Bn 


B w 
m 


0 

0 

B u E h 
m 


in Q , 
on r 1 , 
on r 2 , 


where B m is given by (2.6) (see Figure 1). For simplicity we only con- 
sider the case of axial symmetry, however, all results hold for the general 
three dimensional problem as well as the for the Dirichlet problem. As in 
section 3 we restrict the proofs to the case that the artificial surface 
is the sphere r = r^ We are interested in error estimates on surfaces 
and use the surface norm given by (3.4). We then have 
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Theorem 4.1. 


Given equations (4.2) and m = 1 the error w has the bound 


(4.3) 


w L s < 


'(r x ) 


Kr{> ~ k 


Proof . 

Let w denote the complex conjugate of w. Then from Green's 
theorem and (4.2) we have 


(4.4) 


j' | grad w | 2 dV ~ ^ j |w 2 |dV = J w dA 

n n r„ 


Using the definition of we also have 


(4.5) 


J' w dA = ik ^*|w 2 |dA + h dA 


We then substitute (4.5) into (4.4). Taking the imaginary part of the 
resulting equation yields 


(4.6) w 


(rp -f I "! 2 d A " k /' 


k / wh d A < 


“ “(rp II h H(rp 


Dividing both sides by || w||, . gives the estimate (4.3). 

1 C 

Since h = B^u we know that |h| £ — j where C depends only on g 

r 

and k. Hence, (4.3) is equivalent to 


(4.7) 


w 


' (r l> - kr 2 
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Using the definition of 


we also have 


(4.8) 


in (ri) <» + i>i.ii (ri) + ihi ( , 


£ < 2 + Ill'll, ) £ 4 <2 + H> 

1 r 


In theorem (4.1) we considered the boundary operator B^. We now 
consider the operator B^ but restrict the body T ^ to be a sphere. We 
then have 


Theorem 4.2 

Consider equation (4.2) with m = 2 and F^ the sphere r = 1, 
Then the solution w satisfies the bound 

kr 

(4.9) IM| (r } 1(1 + -^) 

We note that by the construction of B^ we have 

hi <■% , 

r i 

where C depends only on g and k. It follows that 



*»< ri > i 4 

1 r„ 


and (4.9) can be restated as 


(4.10) 





) 


- 24 - 



A similar result holds for the Dirichlet problem. 


Proof 


By separation of variables (assuming axial symmetry) the solution to 
(4.2) has the expansion 


(4 - n) w = £ H (kr)Pp(cos 9) , 

£=0 

where are the spherical Bessel function and P^ are the Legrendre 

polynomials. Hence satisfies the ordinary differential equation 

(4.12) + ~ + (k 2 --^-)H £ = 0 ; A = £(£+l) 

r 

Similarly, by completeness, B m u = h has an expansion 

00 

( 4,;L3 ) h = ^2 h ? P (cos9) 

£=0 

By orthogonality it is sufficient to only consider one term in (4.11). 

We thus have reduced the problem to a one dimensional problem. Let <p = H^(r), 



then 





(4.14a) 


<J>" + | + (k 2 - 

~)<p = 0 , 

r 

• 

(4.14b) 


*’(1) = o 

9 

* 

(4.14c) 

b 2 4> = <t>" 

+ 7 ( f>' + (^- k 2 )(J) - 

2 i k (<f> ' +-^<f>) = h £ 
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We note that for the Dirichlet problem (4.14b) is replaced by 
(f)(1) = 0. For the Neumann problem we define q by q = (f)(1). We then 
introduce a normalized variable 


(4.15) 


^ = 4>/q 


It follows form (4.14) that ip is real. We then take the real and 
imaginary parts of (4.14c) and replace (p by if). Using (4.12) to 
eliminate the second derivatives, we obtain 


ib' + — ib 
r i 



Case I: 



X = 1 ( 1 + 1 ) . 


In this case ip(l) = 1, if) 1 (1) = 0 together with (4.12) imply that 
if), if) 1 , ip” are all positive in a neighborhood of r = 1. This condition 
must persist at least until r >_ r ^ where ip" (r*) < 0, ip f (r^ f ) = 0 
and ip(r^) > 0. However, combining these conditions with (4.12) we see 
that r^ > r^. Since we are only interested in 1 £ r £ r^ we conclude 
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that are positive in the whole interval. In particular ijj(r^) > 0, 


^'(rp > 0. Using (4.16) we have 


2 ~ I £' 

ip(r 1 ) <ip’ + — V 1 I I 1 2kfq’ 


Since 4) = q ^ we have 


(4.17) 


|<Kr 1 )| 


kr 1 h 0 

— I— I 

k 


which is a stronger estimate than (4.9) 
We now consider 


Case II: 


\ < * 2 

r i 


We now solve (4.16) for ip(r^) and obtain 


(4.18) 


iK r i) (rz ~ \ ~ 2r2) 


r l r l 


A 2 

Since £ k we have that 
r l 


i_A_ _2_ 
2 " 2 
r l r l 


2k 2 1 > k 2 + \ 

r l 


Hence, using (4.18) and <j> - q ^ we have 


2,.. 2 , IM , 


k (1+ “0> ~ 

r^ 1 H 


_ < d+^ , 

I ~ I q rk 
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or 


(4.19) 



> 


which again is a stronger estimate than (4.9). We finally point out that 
(4.9) and (4.16a) imply 

(4.20) IU’(r 1 )|| < 

since h = 0 (tt) • 
r i 

The results of Theorems 4.1 and 4.2 can be stated as 


o Kr -i 

i+ ^ (1+ ih =0 ^> 


(4.21) 


W !l(r 1 )’ II W rH(r 1 ) = 0(1/r f 1) 


for m = 1 and 2. We now extend these results from the outer surface 
r = r^ to the body . Introducing the notation 


w|| (r \ = //M 2d A 

v v r 


for the surface norm; we then have 


Theorem 4.2 

For the problem (4.2) we have 


(4.22) 


II w H(r 1 ) - c 



where C depends only on k and 

Note: (1) An entirely analogous theorem holds for 

Dirichlet problem is considered. 


w , v when the 
n'Krp 
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(2) C is continuous in k and the result holds for k = 0. 

(3) Since u = 0(-^) the relative error expressed in (4.22) is 
smaller than that expressed in (4.21). 


Proof 

Let G (p , q) denote the free space Green's function 

e ±k|p-q| 

(4.23) " 4w| P -ql — * 

Using Green's theorem applied to w in the region ft (see figure 1) we 
have for p g ft (making use of (4.2b) 


(4.24) 


w(p) = f G(p,q)w n (q) - w(q)G fl (p,q)dA c 
i Ji a 




+ 


L 


w(q)G 


n 


(p > q) d A q 


9 


where the normal on I* points toward the exterior and the normal on 
| q| = r^ is in the direction of increasing r. Upon letting p approach 
and using the standard jump relations of potential theory ([15]) we obtain 


(4.25) I(w) = f (G(p,q)w n (q) - w(q)G n (p,q))dA £ 

,J . o 


where I is the integral operator (see [9]) 


I(w) 


w(p) 

2 


q 



e r. 


G n (p,q)w(q)d A . 

q 
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It is clear that the right hand side of (4.25) satisfies (4.22). It is 
shown in [9] that I(w) is invertible except for a discrete set of values 
of k and this establishes (4.22) except for these interior resonances. 

At the interior resonances we replace G by the modified Green's function 
of Ursell ([32]) and use the same proof. 


V. Numerical Results 

We consider the three dimensional Helmholtz equation in spherical 
coordinates with axial symmetry. Using r and 0 as coordinates the 
equation becomes 


(5.1) 


3 2 u 2 9u 

dr 2 r 


r 2 sin 0 


9 / j a 

30 (sin 0 30 ) 


+ k u = 0 


The first two radiation boundary conditions are 


(5.2a) 


3 , 0 = - iku + - = 0 

1 9r r 


and 


(5.2b) B,u = + (£- 2 i k) + (| - 4 i k) J - k 2 u = 0 


9r 


9 2 u 

To eliminate — !r, (5.2b) can be rewritten as 
9r 


(5.2c) 


9 9u . 

3u .1 . . 39 ( 3e sln6) 


2r (— - i k )sin 0 


We only consider the case that the body is a sphere, given by r = r^. 

The radiation conditions are specified on the sphere r = r^. On the body 
Neumann data 
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(5.3) 


on 


r 


3u . 
^ = h 


= r r 


is given. 

The Helmholtz equation is solved by a finite element code using linear 
elements. The radiation boundary condition is also written in weak form 
following an integration by parts. Details are given in the appendix. 
Richardson extrapolation in the theta direction is used for all results. 

This substantially increases the accuracy of the computed solutions, 
especially for the higher frequencies. The limiting factor in the method 
is the storage requirements of the banded Gauss solver. The use of higher 
order elements or iterative techniques would reduce the storage requirements. 

We first consider acoustic scattering by a point source at a fixed 
axial point q. We then have 



a) 

Au 

+ k^u = 5 (p-q) 

> 

1-1 

o 

A 

A 

8 

(5.4) 

b) 



5 

II 

O 

v» 


c) 

3u 

9r 

- i k u = o (— ) 
r 

5 

r °° 


Here, p denotes the dependent variable and q is the source point. (5.4) 
implies that the scatterer is hard, i.e., the normal acoustic velocity is zero 
on the body. 

The singularity in (5.4) is eliminated by introducing v given by 

e ik l P-q| 

V = U - 

4tt p-q 
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The equations for v are then 


a) A v + k^v = 0, 



c) - ik v = o(-^) , 


r -> co. 


To solve this system according to the previously described procedure, we 
replace (5.5c) by 


(5.5d) 


B-jV = 0, 


r = r. 


or 


(5.5e) B 2 v = 0, r = Ty 

The "exact 11 solution was generated by using an integral equation code 

with a fine grid. This code had been previously checked by comparisons 

with analytic solutions presented in [6] to verify its accuracy. 

2 

In Table I we present the relative surface L errors for various k. 
The data in this table was obtained with a fixed number of points (N=5) 
in the radial direction. For these computations the body has radius r^ = 
The source is located on the axis at r = .6. Similar results have been 
obtained for a wide range of source positions. 

The results given in Table I indicate that for the problem considered 
the error due to the artificial boundary decreases as k increases. Speci- 
fically the grid in the normal direction can be chosen independent of k. 
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One can also bring in the artificial surface so that it coincides with the 
inner boundary. In this case only an ordinary differential equation needs 
to be solved. However, this resulted in substantial errors and accurate 
solutions could not be obtained. 

The results of Table I show that the first order boundary condition 
(5.2a) is not good enough for many computations. Except for the lowest 
of frequencies 10 percent accurate solutions could not be obtained using 
( 5 . 2 a) because of storage difficulties that arise from using a large r^. 

We next consider the Helmholtz equation with a variable k. We choose 
k in (5.4a) as 


r- r r 




k = < 


r o± r ± r o + * 15 


Tq + .15 £ r. 


We choose r^ = i and consider a sequence of outer boundary positions 
r^. The "exact" solution is generated by choosing the outer boundary 
sufficiently far away so that the solutions obtained by using the first or 
second order boundary conditions differed by less than 3 percent. The solu- 
tion obtained by using at this r^ was taken as the "exact" solution. 

Since k is variable inside the region the Green 1 s function is not 

known. Hence, one can not use Green’s formula to calculate the far field 

9 u 3u 

solution given u and at r = r^. Instead u and on the outer 

boundary can be used to find the far field solution, since k is constant 

3 u 

beyond the outer boundary. The value of 7 ^- on the outer boundary can be 
calculated from u and its tangential derivatives by using the boundary 
conditions (5.2a or 5.2c). The tangential derivatives can be eliminated 
by integration by parts in the Green’s formula. 
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In figure 2 we plot the relative surface errors over the outer 
surface for different values of r^. The improvement achieved by using 
(5,2c) rather than (5.2a) is evident. 

As the final example we consider Laplace's equation. The exact solu- 
tions are generated by imposing Neumann data corresponding to a monopole 
or dipole centered at an axial point p g inside r^ and displaced from 
the origin of coordinates. The solution is 


(5.7a) 





and 


(5.7b) 


u 


cos 8 T 



9 


for the monopole and dipole respectively. Here 0' denotes the polar 
angle in a polar coordinate system centered at p g . 

It is well known that general solutions to Laplace's equation in 
exterior domains can be expressed as a superposition of surface monopoles 
or surface dipoles ([15]). Hence, the model problems (5.7) are relevant 
to realistic problems especially when p g is chosen near the body surface. 

In table 2 we present the surface errors, over the inner surface 
(r = 2 ) 9 f° r several different cases. The solution types M and D denote 
the monopole (5.7a) or dipole (5.7b) solution respectively. 

The first order condition is exact for a monopole centered at the origin 
while the second order boundary condition is exact for both monopoles and 
dipoles centered at the origin. When the second order condition is used for 
the displaced dipole (5.7b) it is more accurate than the first order condition 
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is for a displaced monopole (5.7a). For small p the dipole has no 1/r 
contribution and the first order condition is expected to be inaccurate. 

These results again confirm the substantial improvements that can be ob- 
tained by the use of (5.2c). The decay rate of the error as a function 
of r^, as predicted in Theorem (3.1) has also been computationally verified. 



WKP5QK 



Figure 2.- Errors for variable coefficient problem. 



TABLE 1 


Relative surface L 


errors for the Helmholtz equation 


kr 0 

r l 

m 

error 

11 

.645 

i 

18.6 

11 

.645 

2 

5.4 

23.6 

.575 

1 

31.5 

23.6 

.575 

2 

6.5 


* 
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TABLE 2 


Relative surface L 


errors for the Laplace equation 


Solution type 

p s 

r l 

m 

error 

M 

.2 

.57 

1 

14.7 

M 

.2 

.57 

2 

1.1 

M 

.2 

. 64 

1 

9.4 

M 

.2 

.64 

2 

0.6 

D 

.2 

.57 

i 

65.0 

D 

.2 

.57 

2 

8.9 

D 

.2 

.64 

1 

38.3 

D 

.2 

.64 

2 

4.9 

M 

.4 

.60 

1 

25.4 

M 

.4 

.60 

2 

4.1 

M 

.4 

.65 

1 

17.7 

M 

.4 

.65 

2 

2.7 

D 

.4 

.60 

1 

35.4 

D 

.4 

.60 

2 

10.5 

D 

.4 

.65 

1 

21.6 

D 

.4 

.65 

2 

7.2 
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APPENDIX 


In this section we consider the problem of implementing the boundary 
conditions corresponding to and From (2.6) we see that these 

conditions have the general form 

(A. 1) BjU = + Y-jU = 0, 

and 

2 

(a.2) B 2 u -ri + 6 ilF + v°- 

3r 

Here a^, and y^ are complex functions depending only on k and r 

and not on the angular variables which are given explicitly in (5.2a) and 
(5.2b). 

For simplicity we will consider only axially symmetric problems, so 
that the Helmholtz equation 

(A. 3) A u + k^u = 0, 


may be written as 


(A. 4) 


3 2 u 2 3u L_ 

a 2 + r 3r + 2 . , 

3r r sin i 


3 . 3ik . 2 _ 

0 Q (sin 0 90 ) + k u 0 . 


The presence of the 


3 2 u 

3r 2 


term in (A.2) is nonstandard and difficult to 


implement directly. Since the boundary will be the sphere r = r^ 


it is convenient to eliminate the 


3f» 

3r 2 


term in (A.2) by using (A. 4) 


so that only tangential second derivatives appear in the boundary condition. 
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This can also be done for arbitrary shapes of the artificial boundary. The 
result can be express in the form 


(A. 5) 


V = y 2 I? + li^e !e (sln 0 !e> + a 2 u = °» 


where o^, 3 2 > Y 2 depend only on k and r and not on 0 and are given 
explicitly in (5.2c). The functional form of (A. 5) also includes the 
boundary operator B 1 (with 6 2 = 0). 

The computational problem is to solve (A. 3) in a region fi exterior 
to an inner boundary and interior to the sphere r = r^ On r = ^ 

we impose (A. 5) while on 1^ we may impose either Dirichlet or Neumann 
data. For concreteness we assume that Neumann data is specified, i.e. 

(A ’ 6 ) = g(r,9) on r r 

The numerical method employed is a Galerkin finite element technique which 
we now describe. 

Consider the following weak formulation of the problem (A. 3), (A. 5), 

(A. 6). We seek a function u in a Hilbert space H such that 

m 

(A. 7) B(u,v) = F (v) , 


for all v in H where 
m 


(A. 8) B(u,v) - /yivu -Vv- k 2 uv]r 2 sin0 d0 dr + f [^- || || - ^ uv]r 2 sin 0 d0, 

r\ ** 2 2 


r=r 
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and 


4 


" I s 


F(v) = -J g(r,0)vd s. 


The boundary integral terms in (A. 8) come from using (A, 5) and integrating 

2 

by parts in 0. We note that r sin 0 comes from the three dimensional 

volume element. Our approximation problem is to choose a finite dimensional 

subspace S* 1 of H and then seek a xi 1 e such that 

r m 


(A. 9) 


B(u h ,v h ) = F(v h ), 


for all v^ 1 e S h . 

It is easy to see that if u is a sufficiently smooth solution to (A. 7) 

then u will be a solution to (A. 3), (A, 5) and (A. 6). For m = 1 the 

above weak formulation almost leads to a standard finite element problem, 
h 

once S is chosen to be some finite element space. Indeed, for m = 1, 

$2 = 0 and we may choose = H"^(r), i.e., the Sobolev space of functions 

with one distributional derivative. The only complication is the appearance 
2 

of the term (r sin 0) in the integrals. For m = 2, we need more smooth- 
ness on the boundary r = r^ due to the first term in the boundary integral 

in (A. 8). The choice x H^(r=r^) suffices. For either m = 1 

h 

or m = 2 we may choose 5 to be the finite element space defined by 

subdividing into triangular elements and then restricting the function in 

1 ^ 

S to be continuous in Q and linear in each triangle, i.e. a piecewise 

linear finite element space. Then we compute u in the standard manner, 

h h 

i.e. we choose a (local) basis for S , expand u in terms of this basis, 
and let v in (A. 9) range over the basis to obtain a system of algebraic 
equations for the nodal values of u . 
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Numerical experiments confirm that when both the 0 and r grids 

2 

are refined that h convergence rates are obtained. These convergence 
rates were verified for both the first order B^u = 0 and second order 
B^u = 0 boundary conditions. Hence, the second order derivatives that 
appear in did not cause any deterioration in the rate of convergence. 

When only the 0 mesh was refined and the errors were measured on the 
inner boundary a convergence rate of order h was observed. 

The computations presented in this paper were intended to exhibit the 
improvements due to increasing the order of the boundary operators. For 
this reason, it was desired to minimize discretization errors due to the 
finite element approximation. It has been verified that the higher order 
conditions did not affect the rate of convergence of the discretization 


scheme. 
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